home *** CD-ROM | disk | FTP | other *** search
/ Sprite 1984 - 1993 / Sprite 1984 - 1993.iso / src / lib / m / lgamma.c < prev    next >
C/C++ Source or Header  |  1988-07-11  |  3KB  |  148 lines
  1. /*
  2.  * Copyright (c) 1985 Regents of the University of California.
  3.  * All rights reserved.  The Berkeley software License Agreement
  4.  * specifies the terms and conditions for redistribution.
  5.  */
  6.  
  7. #ifndef lint
  8. static char sccsid[] = "@(#)lgamma.c    5.2 (Berkeley) 4/29/88";
  9. #endif /* not lint */
  10.  
  11. /*
  12.     C program for floating point log Gamma function
  13.  
  14.     lgamma(x) computes the log of the absolute
  15.     value of the Gamma function.
  16.     The sign of the Gamma function is returned in the
  17.     external quantity signgam.
  18.  
  19.     The coefficients for expansion around zero
  20.     are #5243 from Hart & Cheney; for expansion
  21.     around infinity they are #5404.
  22.  
  23.     Calls log, floor and sin.
  24. */
  25.  
  26. #include <math.h>
  27. #if defined(vax)||defined(tahoe)
  28. #include <errno.h>
  29. #endif    /* defined(vax)||defined(tahoe) */
  30. int    signgam = 0;
  31. static double goobie    = 0.9189385332046727417803297;    /* log(2*pi)/2 */
  32. static double pi    = 3.1415926535897932384626434;
  33.  
  34. #define M 6
  35. #define N 8
  36. static double p1[] = {
  37.     0.83333333333333101837e-1,
  38.     -.277777777735865004e-2,
  39.     0.793650576493454e-3,
  40.     -.5951896861197e-3,
  41.     0.83645878922e-3,
  42.     -.1633436431e-2,
  43. };
  44. static double p2[] = {
  45.     -.42353689509744089647e5,
  46.     -.20886861789269887364e5,
  47.     -.87627102978521489560e4,
  48.     -.20085274013072791214e4,
  49.     -.43933044406002567613e3,
  50.     -.50108693752970953015e2,
  51.     -.67449507245925289918e1,
  52.     0.0,
  53. };
  54. static double q2[] = {
  55.     -.42353689509744090010e5,
  56.     -.29803853309256649932e4,
  57.     0.99403074150827709015e4,
  58.     -.15286072737795220248e4,
  59.     -.49902852662143904834e3,
  60.     0.18949823415702801641e3,
  61.     -.23081551524580124562e2,
  62.     0.10000000000000000000e1,
  63. };
  64.  
  65. double
  66. lgamma(arg)
  67. double arg;
  68. {
  69.     double log(), pos(), neg(), asym();
  70.  
  71.     signgam = 1.;
  72.     if(arg <= 0.) return(neg(arg));
  73.     if(arg > 8.) return(asym(arg));
  74.     return(log(pos(arg)));
  75. }
  76.  
  77. static double
  78. asym(arg)
  79. double arg;
  80. {
  81.     double log();
  82.     double n, argsq;
  83.     int i;
  84.  
  85.     argsq = 1./(arg*arg);
  86.     for(n=0,i=M-1; i>=0; i--){
  87.         n = n*argsq + p1[i];
  88.     }
  89.     return((arg-.5)*log(arg) - arg + goobie + n/arg);
  90. }
  91.  
  92. static double
  93. neg(arg)
  94. double arg;
  95. {
  96.     double t;
  97.     double log(), sin(), floor(), pos();
  98.  
  99.     arg = -arg;
  100.      /*
  101.       * to see if arg were a true integer, the old code used the
  102.       * mathematically correct observation:
  103.       * sin(n*pi) = 0 <=> n is an integer.
  104.       * but in finite precision arithmetic, sin(n*PI) will NEVER
  105.       * be zero simply because n*PI is a rational number.  hence
  106.       *    it failed to work with our newer, more accurate sin()
  107.       * which uses true pi to do the argument reduction...
  108.       *    temp = sin(pi*arg);
  109.       */
  110.     t = floor(arg);
  111.     if (arg - t  > 0.5e0)
  112.         t += 1.e0;                /* t := integer nearest arg */
  113. #if defined(vax)||defined(tahoe)
  114.     if (arg == t) {
  115.         extern double infnan();
  116.         return(infnan(ERANGE));        /* +INF */
  117.     }
  118. #endif    /* defined(vax)||defined(tahoe) */
  119.     signgam = (int) (t - 2*floor(t/2));    /* signgam =  1 if t was odd, */
  120.                         /*            0 if t was even */
  121.     signgam = signgam - 1 + signgam;    /* signgam =  1 if t was odd, */
  122.                         /*           -1 if t was even */
  123.     t = arg - t;                /*  -0.5 <= t <= 0.5 */
  124.     if (t < 0.e0) {
  125.         t = -t;
  126.         signgam = -signgam;
  127.     }
  128.     return(-log(arg*pos(arg)*sin(pi*t)/pi));
  129. }
  130.  
  131. static double
  132. pos(arg)
  133. double arg;
  134. {
  135.     double n, d, s;
  136.     register i;
  137.  
  138.     if(arg < 2.) return(pos(arg+1.)/arg);
  139.     if(arg > 3.) return((arg-1.)*pos(arg-1.));
  140.  
  141.     s = arg - 2.;
  142.     for(n=0,d=0,i=N-1; i>=0; i--){
  143.         n = n*s + p2[i];
  144.         d = d*s + q2[i];
  145.     }
  146.     return(n/d);
  147. }
  148.